Libraries:
# libraries -------------------
library(ggplot2)
library(vegan)
library(ggvegan)
library(tidyverse)
library(ggordiplots)
Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')
Functions:
# functions -------------------
# 3 dimensions
mali_3d <- metaMDS(rel_mali5,
distance = 'bray',
k = 3,
trymax = 20,
autotransform = FALSE)
## Run 0 stress 0.1272962
## Run 1 stress 0.1275737
## ... Procrustes: rmse 0.03524096 max resid 0.1289734
## Run 2 stress 0.1275741
## ... Procrustes: rmse 0.03471258 max resid 0.1272658
## Run 3 stress 0.127296
## ... New best solution
## ... Procrustes: rmse 0.0001272258 max resid 0.000417742
## ... Similar to previous best
## Run 4 stress 0.1582628
## Run 5 stress 0.1272958
## ... New best solution
## ... Procrustes: rmse 0.000587606 max resid 0.001999951
## ... Similar to previous best
## Run 6 stress 0.1275741
## ... Procrustes: rmse 0.03499637 max resid 0.1277436
## Run 7 stress 0.1272965
## ... Procrustes: rmse 0.0008573844 max resid 0.002921016
## ... Similar to previous best
## Run 8 stress 0.1275736
## ... Procrustes: rmse 0.03451709 max resid 0.1262122
## Run 9 stress 0.1275735
## ... Procrustes: rmse 0.03479077 max resid 0.127088
## Run 10 stress 0.1272959
## ... Procrustes: rmse 0.0004557986 max resid 0.001556358
## ... Similar to previous best
## Run 11 stress 0.1272965
## ... Procrustes: rmse 0.0008645533 max resid 0.002945962
## ... Similar to previous best
## Run 12 stress 0.1275742
## ... Procrustes: rmse 0.03502329 max resid 0.1277945
## Run 13 stress 0.1272961
## ... Procrustes: rmse 0.0002420751 max resid 0.0007293543
## ... Similar to previous best
## Run 14 stress 0.1275734
## ... Procrustes: rmse 0.03460257 max resid 0.1264966
## Run 15 stress 0.127574
## ... Procrustes: rmse 0.03434109 max resid 0.1256504
## Run 16 stress 0.1272962
## ... Procrustes: rmse 0.0007275598 max resid 0.002468901
## ... Similar to previous best
## Run 17 stress 0.1275736
## ... Procrustes: rmse 0.03484616 max resid 0.1272575
## Run 18 stress 0.1272959
## ... Procrustes: rmse 8.512657e-05 max resid 0.000239125
## ... Similar to previous best
## Run 19 stress 0.1272963
## ... Procrustes: rmse 0.0007570094 max resid 0.002579365
## ... Similar to previous best
## Run 20 stress 0.1275734
## ... Procrustes: rmse 0.03459637 max resid 0.1264763
## *** Best solution repeated 8 times
# converges, stress of 0.127 - 0.128
set.seed(2)
# 2 dimensions
mali_2d <- metaMDS(rel_mali5,
distance = 'bray',
k = 2,
try = 40,
trymax = 40,
autotransform = FALSE)
## Run 0 stress 0.2013098
## Run 1 stress 0.2421891
## Run 2 stress 0.2164813
## Run 3 stress 0.2248403
## Run 4 stress 0.2013098
## ... New best solution
## ... Procrustes: rmse 7.336217e-05 max resid 0.0002404534
## ... Similar to previous best
## Run 5 stress 0.2166909
## Run 6 stress 0.2013098
## ... Procrustes: rmse 2.20954e-05 max resid 7.187511e-05
## ... Similar to previous best
## Run 7 stress 0.2299679
## Run 8 stress 0.2334747
## Run 9 stress 0.2067192
## Run 10 stress 0.2063475
## Run 11 stress 0.2248403
## Run 12 stress 0.2052525
## Run 13 stress 0.2013098
## ... Procrustes: rmse 8.623799e-05 max resid 0.0002840637
## ... Similar to previous best
## Run 14 stress 0.2248403
## Run 15 stress 0.2013098
## ... Procrustes: rmse 7.687422e-05 max resid 0.0002349284
## ... Similar to previous best
## Run 16 stress 0.2385195
## Run 17 stress 0.2378598
## Run 18 stress 0.2013098
## ... Procrustes: rmse 7.160904e-05 max resid 0.0002354649
## ... Similar to previous best
## Run 19 stress 0.240332
## Run 20 stress 0.271914
## Run 21 stress 0.2063472
## Run 22 stress 0.2013098
## ... Procrustes: rmse 3.826719e-05 max resid 9.76465e-05
## ... Similar to previous best
## Run 23 stress 0.2337596
## Run 24 stress 0.2518641
## Run 25 stress 0.2401184
## Run 26 stress 0.2013098
## ... Procrustes: rmse 6.309312e-06 max resid 1.438008e-05
## ... Similar to previous best
## Run 27 stress 0.2431912
## Run 28 stress 0.2013098
## ... Procrustes: rmse 2.79146e-06 max resid 8.73469e-06
## ... Similar to previous best
## Run 29 stress 0.2076391
## Run 30 stress 0.2013099
## ... Procrustes: rmse 0.0001601721 max resid 0.0005288505
## ... Similar to previous best
## Run 31 stress 0.2013098
## ... Procrustes: rmse 2.72951e-05 max resid 6.009303e-05
## ... Similar to previous best
## Run 32 stress 0.2185542
## Run 33 stress 0.2248403
## Run 34 stress 0.2248403
## Run 35 stress 0.239287
## Run 36 stress 0.2248403
## Run 37 stress 0.2076443
## Run 38 stress 0.2394146
## Run 39 stress 0.2013098
## ... Procrustes: rmse 1.065455e-05 max resid 3.451836e-05
## ... Similar to previous best
## Run 40 stress 0.2013098
## ... Procrustes: rmse 2.107096e-05 max resid 6.965666e-05
## ... Similar to previous best
## *** Best solution repeated 12 times
# stress still too high to represent in 2 dimensions, lowest is ~ 0.2
rm(mali_2d)
ma.veg_enr <- envfit(mali_3d, ma.meta.data_veg)
ma.veg_spc.envr <- envfit(mali_3d, rel_mali5)
# site data -------------------
ma.site.scrs <- as.data.frame(scores(mali_3d, display = "sites"))
ma.site.scrs <- cbind(ma.site.scrs, Treat_Type = ma.meta.data_veg$Treat_Type)
ma.site.scrs <- cbind(ma.site.scrs, Region = ma.meta.data_veg$Region)
# species data -------------------
ma.spp.scrs <- as.data.frame(scores(ma.veg_spc.envr, display = "vectors"))
ma.spp.scrs <- cbind(ma.spp.scrs, Species = rownames(ma.spp.scrs))
ma.spp.scrs <- cbind(ma.spp.scrs, pval = ma.veg_spc.envr$vectors$pvals)
sig.ma_spp.scrs <- subset(ma.spp.scrs, pval<=0.05)
print(sig.ma_spp.scrs)
## NMDS1 NMDS2 Species pval
## GAPR -0.43281215 -0.36593981 GAPR 0.028
## VAAN -0.58546249 0.23635680 VAAN 0.003
## QUCO 0.47196034 0.31546257 QUCO 0.009
## PIRI -0.09465907 0.63991646 PIRI 0.001
## KAAN -0.64342045 0.04947569 KAAN 0.001
## MELI -0.21624713 -0.54654802 MELI 0.003
## RUSP -0.62454425 0.06726456 RUSP 0.004
## SMGL 0.06314658 0.74452749 SMGL 0.001
## QUPR -0.58228890 -0.08839473 QUPR 0.012
# create species vectors for indicator species -------------------
mali.ind.spp <- ma.spp.scrs %>%
filter(Species %in% c('QUPR', "KAAN", "RUSP", "SMGL"))
# envr data -------------------
ma.envr.scrs <- as.data.frame(scores(ma.veg_enr, display = "vectors"))
ma.envr.scrs <- cbind(ma.envr.scrs, env.variables = rownames(ma.envr.scrs))
ma.envr.scrs <- cbind(ma.envr.scrs, pval = ma.veg_enr$vectors$pvals)
ma.sig_envr.scrs <- subset(ma.envr.scrs, pval<=0.05)
print(ma.sig_envr.scrs)
## NMDS1 NMDS2 env.variables pval
## Moss -0.1485843 0.5469886 Moss 0.011
## BA_HA 0.5191909 -0.3902654 BA_HA 0.004
## PIRI.BA_HA 0.5939896 -0.4423759 PIRI.BA_HA 0.001
ma.sig_envr.scrs <- ma.sig_envr.scrs %>%
filter(env.variables != "BA_HA")
I’m only going to keep one BA measure; PIRI has a lower p value, so thinking that one, but open to other ideas.
Using the first two axes, which represent the most variation in the ordination space
library(plotly)
fort.mali3 <- fortify(mali_3d) %>%
filter(score == "sites")
plot_ly(x = fort.mali3$NMDS1, y = fort.mali3$NMDS2, z= fort.mali3$NMDS3, type = "scatter3d", mode = "markers", color = ma.meta.data_veg$Treat_Type, )
# Here is a 3D plot, just to show
stressplot(mali_3d)
ordicloud(mali_3d)
# based on the idea that NMDS scores express the highest variation in first axis and on down, I will plot the 1st two and take a look -------------------
# control = yellow, fallrx = dark green, harvest = beige, mowrx = light green, spring rx = red
colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")
# create ellipse plot -------------------
mali.plot.e <- gg_ordiplot(mali_3d, groups = ma.meta.data_veg$Treat_Type)
# pull out ellipse data -------------------
ellipse.mali <- fortify(mali.plot.e$df_ellipse) %>%
group_by(Group)
# base plot -------------------
ma.nmds.plot <- ggplot(ma.site.scrs, aes(x=NMDS1, y=NMDS2))+
geom_point(aes(NMDS1, NMDS2, colour = factor(ma.site.scrs$Treat_Type), shape = factor(ma.site.scrs$Region)), size = 2)+
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
coord_fixed()+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type", shape = "Region")+
theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10))
# plot with ellipse -------------------
ma.nmds.plot +
geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group), linejoin = "round", show.legend = FALSE)+
labs(title = "Coastal barrens ordination \n (Axes 1 & 2)")+
theme(plot.title=element_text(hjust=0.5, size=20))
# add significant species -------------------
ma.nmds.plot +
geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group), linejoin = "round", show.legend = FALSE)+
geom_segment(data = sig.ma_spp.scrs, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", lwd = 0.3)+
ggrepel::geom_text_repel(data = sig.ma_spp.scrs, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
labs(title = "Coastal barrens ordination with \n significant species (Axes 1 & 2)")+
theme(plot.title=element_text(hjust=0.5, size=20))
# add indicators species -------------------
ind.sp.mali_color <- c("QUPR" = "#81A88D", "KAAN" = "#81A88D", "SMGL" = "#A2A475", "RUSP" = "#81A88D")
ma.nmds.plot +
geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group), linejoin = "round", show.legend = FALSE)+
geom_segment(data = mali.ind.spp, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = ind.sp.mali_color, lwd = 0.75)+
ggrepel::geom_text_repel(data = mali.ind.spp, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
labs(title = "Coastal barrens ordination with \n indicator species (Axes 1 & 2)")+
theme(plot.title=element_text(hjust=0.5, size=20))
# add significant envr factors -------------------
ma.nmds.plot +
geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group), linejoin = "round", show.legend = FALSE)+
geom_segment(data = ma.sig_envr.scrs, aes(x = 0, xend = NMDS1, y = 0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")))+
ggrepel::geom_text_repel(data = ma.sig_envr.scrs, aes(x=NMDS1, y=NMDS2, label = env.variables), cex = 4, direction = "both", segment.size = 0.25) +
labs(title = "Coastal barrens ordination with \n environmental vectors (Axes 1 & 2)")+
theme(plot.title=element_text(hjust=0.5, size=20))
ma.nmds.plot +
geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group), linejoin = "round", show.legend = FALSE)+
geom_segment(data = sig.ma_spp.scrs, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "grey10", lwd = 0.3)+
ggrepel::geom_text_repel(data = sig.ma_spp.scrs, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
geom_segment(data = mali.ind.spp, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = ind.sp.mali_color, lwd = 0.75)+
ggrepel::geom_text_repel(data = mali.ind.spp, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
labs(title = "Coastal barrens ordination with significant \n species & environmental factors (Axes 1 & 2)")+
theme(plot.title=element_text(hjust=0.5, size=20))
ma.nmds.plot +
geom_path(data = ellipse.mali, aes(x=x, y=y, color=Group), linejoin = "round", show.legend = FALSE)+
geom_segment(data = mali.ind.spp, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = ind.sp.mali_color, lwd = 0.75)+
ggrepel::geom_text_repel(data = mali.ind.spp, aes(x=NMDS1, y=NMDS2, label = Species), cex = 3, direction = "both", segment.size = 0.25)+
geom_segment(data = ma.sig_envr.scrs, aes(x = 0, xend = NMDS1, y = 0, yend=NMDS2), arrow = arrow(length = unit(0.25, "cm")))+
ggrepel::geom_text_repel(data = ma.sig_envr.scrs, aes(x=NMDS1, y=NMDS2, label = env.variables), cex = 4, direction = "both", segment.size = 0.25) +
labs(title = "Coastal barrens ordination with indicator \n species & environmental factors (Axes 1 & 2)")+
theme(plot.title=element_text(hjust=0.5, size=20))
summary(as.factor(ma.meta.data_veg$Treat_Type))
## Control FallRx Harvest MowRx SpringRx
## 5 3 6 3 6
ma.tt <- adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
print(ma.tt)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 4 1.9569 0.31903 2.1082 0.001 ***
## Residual 18 4.1769 0.68097
## Total 22 6.1339 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(pairwiseAdonis)
pairwise.adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
## $parent_call
## [1] "rel_mali5 ~ Treat_Type , strata = Null , permutations 999"
##
## $FallRx_vs_SpringRx
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.21549 0.10956 0.8613 0.614
## Residual 7 1.75138 0.89044
## Total 8 1.96687 1.00000
##
## $FallRx_vs_MowRx
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.48859 0.46035 3.4122 0.1
## Residual 4 0.57276 0.53965
## Total 5 1.06134 1.00000
##
## $FallRx_vs_Harvest
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.40379 0.20243 1.7767 0.109
## Residual 7 1.59094 0.79757
## Total 8 1.99473 1.00000
##
## $FallRx_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.39033 0.23712 1.8649 0.12
## Residual 6 1.25581 0.76288
## Total 7 1.64614 1.00000
##
## $SpringRx_vs_MowRx
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.65697 0.28336 2.7678 0.015 *
## Residual 7 1.66151 0.71664
## Total 8 2.31848 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $SpringRx_vs_Harvest
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.4952 0.15598 1.8481 0.061 .
## Residual 10 2.6797 0.84402
## Total 11 3.1749 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $SpringRx_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.16729 0.0666 0.6422 0.778
## Residual 9 2.34457 0.9334
## Total 10 2.51186 1.0000
##
## $MowRx_vs_Harvest
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.72529 0.32578 3.3823 0.014 *
## Residual 7 1.50107 0.67422
## Total 8 2.22636 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $MowRx_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.84147 0.41918 4.3303 0.012 *
## Residual 6 1.16594 0.58082
## Total 7 2.00741 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Harvest_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.58418 0.21102 2.4072 0.023 *
## Residual 9 2.18412 0.78898
## Total 10 2.76830 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"
#significant: SpringRx vs. MowRx; SpringRx vs. Harvest (borderline at 0.051); MowRx vs. Harvest; MowRx vs Control; Harvest vs Control
library(labdsv)
ma.meta.data_veg$Treat_Group <- NA
ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "FallRx", 1, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "MowRx", 2, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "SpringRx", 3, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Harvest", 4, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Control", 5, ma.meta.data_veg$Treat_Group)
mali_isa <- indval(mali5_avg, ma.meta.data_veg$Treat_Group)
summary(mali_isa)
## cluster indicator_value probability
## QUPR 2 0.8549 0.001
## KAAN 2 0.7534 0.012
## RUSP 2 0.7097 0.006
## SMGL 4 0.6896 0.009
##
## Sum of probabilities = 4.689
##
## Sum of Indicator Values = 6.68
##
## Sum of Significant Indicator Values = 3.01
##
## Number of Significant Indicators = 4
##
## Significant Indicator Distribution
##
## 2 4
## 3 1
gr <- mali_isa$maxcls[mali_isa$pval <= 0.05]
iv <- mali_isa$indcls[mali_isa$pval <= 0.05]
pv <- mali_isa$pval[mali_isa$pval <= 0.05]
fr <- apply(mali5_avg > 0, 2, sum)[mali_isa$pval <= 0.05]
fridg <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
fridg <- fridg[order(fridg$group, -fridg$indval),]
print(fridg)
## group indval pvalue freq
## QUPR 2 0.8549332 0.001 6
## KAAN 2 0.7534400 0.012 10
## RUSP 2 0.7096551 0.006 7
## SMGL 4 0.6896051 0.009 6
# 1 is FallRx
# 2 is MowRx
# 3 is SpringRx
# 4 is Harvest
# 5 is Control